Anisotropic surface tension of buckled fluid membrane 
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Thin solid sheets and fluid membranes exhibit buckling under lateral compression. Here, it is 
revealed that buckled fluid membranes have anisotropic mechanical surface tension contrary to solid 
sheets. Surprisingly, the surface tension perpendicular to the buckling direction shows stronger 
dependence on the projected area than that parallel to it. Our theoretical predictions are supported 
by numerical simulations of a meshless membrane model. This anisotropic tension can be used 
to measure the membrane bending rigidity. It is also found phase synchronization occurs between 
multilayered buckled membranes. 
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I. INTRODUCTION 

Buckling and crumpling of thin solid sheets and strips 
1,0| are commonly seen in our daily life (e.g., for fruits 
3|, paper [1, and polyster film 0) as well as on 
the molecular scale (e.g., for viruses [3] and atomistic 
graphene sheets The shapes produced by the buck- 
ling of elastic sheets have fascinated many physicists. In 
1691, J. Bernoulli proposed the problem of a simple bent 
beam or rod "elastica" , which led to the development of 
the calculus of variation and the elliptic function Q . The 
curve on a plane for minimum bending energy J c^ds 
with constant total length L is expressed by elliptic func- 
tions, where c is the curvature and s is the arc length. 
The theory of elastica was recently extended to twisted 
rods to describe the shape of supercoiled DNAs [Toi] and 
it has been em plo yed to draw smooth surfaces in com- 
puter graphics pl| . 

The buckling has been intensively investigated for 
Langmuir monolayers on air- water interface ,12 17]. The 
buckling develops to collapse or fold of the monolayers 
into water. For a fluid membrane, the balance between 
gravity and membrane bending energy determines the 
buckling wavelength [l^. Recently, buckling transition 
was also observed in a fluid bilayer membrane in simu- 
lations (l8l - [20| . For the bilayer membrane, the effects of 
gravitation and membrane dissolution to the supporting 
water are negligible. Thus, the bilayer membrane is a 
simpler system so it is suitable to study the buckling of 
fluid membranes in details. 

Recently, the stress of torque tensors in fluid mem- 
branes were derived [111, (2^. When the membrane is 
curved, the mechanical surface tension is anisotropic and 
deviated from the thermodynamic surface tension (en- 
ergy to create a unit membrane area). For a tubular 
membrane, the axial stress is finite while the azimuth 
stress is zero. 

In this paper, we report on the shape and surface ten- 



sion of buckled fluid membranes using an analytical the- 
ory and numerical simulation. Since the shape of buckled 
membrane can be analytically derived, the anisotropy of 
the mechanical tensions can be investigated in details. 
Buckling is one of the triggers for breaking membranes. 
For example, under shear flow, the formation of multi- 
23l . 2^ is considered to be induced by 
Shear suppression of the thcr- 



lamellar vesicles 
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buckling instability [2£ 
mal undulations of the membrane and the resulting re- 
duction of the excess area induces buckling. Our study 
revealed that buckling produces large anisotropy in the 
mechanical surface tension. This anisotropy may play a 
role in the membrane stability. 



II. SIMULATION METHOD 

We employ a solvent-free meshless membrane model 
(20I [26j to simulate buckling. The fluid membrane is rep- 
resented by a self-assembled one-layer sheet of particles. 
The particles interact with each other via the potential 
U — e(C/icp + C^att) + kaUa, which consists of a repulsive 
soft-core potential U^op with a diameter cr, an attractive 
potential C/att, and a curvature potential C/q. The details 
of the model are given in Appendix. 

We simulate single- and multi-layer fluid membranes 
by Brownian dynamics in the NVT ensemble with peri- 
odic boundary conditions in a rectangular box with side 
lengths Lx, Ly, and L^. The buckling is chosen in the 
X direction. The buckling occurs in the longest wave 
length Lx in the x direction. When is gradually re- 
duced, the buckled membrane with large amplitude keep 
its direction along x axis even at < Ly. The bend- 
ing rigidity Kcv of the membrane shows a linear depen- 
dence on the parameter ka- We calculate Kcv from the 
thermal undulations of the planar membrane [l^j . Three 
typical bending rigidities are chosen for the simulations: 
Kc^/ksT = 9 ± 0.2, 21 ± 0.5, and 44 ± 1 for kc,/kBT = 5, 
10, and 20, respectively, where k^T is the thermal energy. 
The surface tensions are calculated using "/^ = —PxxLz 
and 7j, = —PyyLz, with the diagonal components of 
the pressure tensor Paa = (NksT — X^i*^*^^)/^ for 
a e {x,y, z}, since Pzz — for solvent-free models. The 
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numerical errors are estimated from three or more inde- 
pendent runs. 



III. SINGLE BUCKLED MEMBRANE 

First, we consider the buckhng of an isolated planar 
membrane. Figure [T] shows buckled membranes at small 
projected areas A^y — L^Ly. When the thermal undula- 
tions are neglected, the shape of the membrane is given 
by the energy minimum of the bending energy f cv of the 
membrane with area constraints. For the buckled mem- 
brane, it is given by 



"io 2 \ds) 



ds 



(1) 



LLy, where s and 9 are 



with constant intrinsic area A 
the arc length and the tangential angle on xz plane, re- 
spectively. When we consider the membrane compressed 
by a constant force XLy in the x direction, it is an elas- 
tica problem to minimize Ftotai = Fc-v + XLxLy. Also, 
the force XLy can be considered as a Lagrange multiplier 
to fix the length L^-. Euler's equation gives the shape 
equation 



cos(6') 



J, 



(2) 



where the characteristic length b — 



/A and ^max is 
Then, the arc 



the maximum tangential angle 0, 
length s is written as s = 6 F{tf, k), with k = sin(0niax/2) 
and sin(6'/2) = ksm{(p), where F{(p, k) is the elliptic inte- 
CTal of the first kind with the elliptic modulus < < 1 
[291] . The modulus k is determined by the total arc length 
L from 



(3) 



where K{k) = F{ti/2^ k) is the complete elliptic integral 
of the first kind. The shape of the buckled membrane is 
expressed by 



X ^ 2b E{am{s/b,k),k) 
z — 2kb cn(s/6, k), 



(4) 
(5) 



where E{a^k), am(a. A;), and cn(a, fe) are the elliptic in- 
tegral of the second kind, Jacobi amplitude, and Jacobi 
elliptic function, respectively [29]. Equations ^ and ([5]) 
reproduce the buckled shape of the simulation very well 
(see Fig. [T]). The thermal fluctuations give small un- 
dulations of the simulated membrane around the energy 
minimum shape (solid curve). Strongly buckled mem- 
brane with 6'max > 7''/2 (fc^ > 1/2) has a $7 shape as 
shown in Fig. [ijc). It is called class 4 of Euler's elastica 
i. 

Next, we derive the surface tension for fixed and Ly 
using elliptic functions. We only consider the mechanical 
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FIG. 1: (Color online) Sliced shape of buckled fluid mem- 
branes at Lx = 35cr and Kcv = 21k^T. (a) A^y/A = 0.75, 
Ly = 24cr, and N = 800. (c) A^y/A = 0.375, Ly = 48cr, and 
= 3200. The solid curves are given by Eqs. @ and (O. 
The closed symbols (•, A,o) represent membrane particles in 
the simulation for three sequential sliced snapshots with time 
interval 2000r. The curve and points are shifted in the x di- 
rection to overlap at {x,z) = (0.251/3,, 0). (b) The 3D image 
of the simulated membrane shown with closed circles in (a). 



surface tension here. The modulus k is determined from 
the length ratio. 



Lx 
~L 



2E(k) 
K[k) 



1. 



(6) 



The stress A is a variable given by A = \^K,f^^(K(}i) j F)^ ■ 
The bending energy of the buckled membrane is given by 

Fev = A{(2fc2 - l)L-f L,} 

= ^^^^[{k^ -X)K{kf -E(k)K(k)\ (7) 

The surface tensions in the x and y directions are given 
by 



dF,, 



lx 

Iv 



Ly dLx 



= -A, 



LxdLy 



-X 



2F 

^± cv 

LxLy 



(8) 
(9) 



respectively. The surface tension ^x is balanced by the 
compressed stress A as expected. However, the surface 
tension 7j, has the additional term 2Fcv/ LxLy. This 
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FIG. 2: (Color online) Area Axy dependence of the surface 
tension for Kcv = 21fcBT and Ly = 24(t (A'^ = 800) or Ly = 
48(7 (N = 3200). The symbols (o, A) and (□, V) represent 
the surface tension 7^ and yy in the simulation, respectively. 
The solid lines are given by Eqs. ([SJ and ©. 



anisotropy can be also derived from the stress tensor de- 
rived in Ref. [2l|, In contrast, solid sheets show 
much weaker correlation between the x and y directions 
because of the shear elasticity. 

The buckling transition points are obtained from the 
condition required to satisfy Eq. ([3]). Since K(k) > 
K{0) = 7r/2, it is written as A > Aq = {2tt/L)^Kc^ 
Therefore, the planar membrane becomes unstable at 
7 = —(2tt/LxYhcv This result is in agreement with that 
estimated from the instability of the lowest Fourier mode 
of the thermal undulations [l^ [3 ■ This coincidence 
is not surprising because the elliptic function reduces to 
a trigonometric function at fc = 0. 

Figure [2] shows the surface tension dependence on Axy 
for constant Ly. At Axy < Axy = lANa"^, the membrane 
is buckled, and then the two surface tensions jx and 
7j, show different values. As Axy decreases, it is found 
that 7j, increases, while jx shows a gradual decrease. In- 
terestingly, the area decrease generates reduction of the 
counter stress in the y direction. It is surprising that 
becomes positive at Axy/Na^ < 0.9 so that the buckled 
membrane prefers to shrink in the y direction. The av- 
erage surface tension 7av = (7^ + 7j/)/2 also increases. 
These simulation results are in good agreement with our 
theoretical prediction given by Eqs. (|S]) and © with the 
area A — Axy The larger membrane starts buckling at 
smaller A = —7, and it has smaller Axy dependence of 'jx 
and 7j,, since Aq oc (compare data for N = 800 and 
3200 in Fig. H]). 

When the aspect ratio Ly/Lx is fixed, the surface ten- 
sions show a different type of Axy dependence (see Fig. 
El). With decreas ing Axy, "fx gradually increases, in con- 
trast to its behavior for constant Ly, and the tension 
difference 7^ — "fx increases weakly. These effects are due 



to an increase in the arc length L = A/Ly for the fixed 
aspect ratio. Thus, the surface tensions are dependent 
on the projected area as well as the aspect ratio. 

When the aspect ratio Ly/Lx is allowed to change 
freely for a fixed projected area Axy = LxLy, the mem- 
brane elongates in the x direction {Lx — )■ 00) in order to 
reduce the membrane bending energy {Ly/L — with 
constant k in Eq. ([7]) ): i.e. Mechanically, the membrane 
pushes the wall in the buckled direction more than that 
in the other direction. This is a qualitative explanation 
of the anisotropic surface tension. Thus, the buckling 
gives the effective shear elasticity to the membrane. 

For small buckling amplitude at 1 — Axy/A'xy ^ 1, the 
surface tensions can be expressed by polynomials. Eq. 
dS]) can be expanded to Lx/L = 1 - fc^ - fc''/8 -f o[k^) for 
small k. From this relation and the expansions of Eqs. 
([8]), and (H)), the surface tensions are expressed as 

TT^Kcv/. Lx^^ , 

7. = '^(5--j, (10) 

The area dependence can be clearly captured by these 
equations. The surface tensions are determined by three 
quantities, Lx/L = Axy/A'^y, A'^y, and These equa- 
tions give a good approximation for Lx/L > 0.6 (compare 
solid and dashed lines in Fig. 

Anisotropic surface tensions are also generated in a 
tubular membrane [sot - ls^ ]. The surface tension in the 
axial direction is given by 7ax = i^cv/ , where R is the 
radius of the membrane tube; while the surface tension 
is zero in the azimuth direction (The average surface ten- 
sion 7av = 7ax/2). The bending rigidity Kcv of the mem- 
brane has been measured usin g th is axial tension in ex- 
periments and simulations j3l|. Similarly, be 
measured from the surface tension of the buckled mem- 
brane. We fit the curves in Fig. Eljb) by Eq. ([TT|) with the 
fit parameters Kcv and A'xy for 0.05 < [A'xy— Axy) / N a'^ < 
0.3. This gives K^y/ksT = 9.4±0.3 (9.3±0.3), 22.3±0.1 
(22.9 ± 0.2), and 47.4 ± 0.1 (48.4 ± 0.2) for constant Ly 
(constant ratio Ly/Lx) at ka/k-oT = 5, 10, and 20, re- 
spectively. These values are in reasonable agreement with 
the bending rigidities estimated from the thermal undu- 
lations of the planar membranes. Compared to the tubu- 
lar membrane, this simulation method is easy to apply to 
explicit solvent systems and bilayer membranes with low 
fiip-flop frequency. In the buckling, the area difference 
between the upper and the lower leaflets of bilayers are 
not changed, since the membrane deformation is sym- 
metric. The solvent is not enclosed by the membrane, so 
the solvent volume is conserved when the volume of the 
simulation box is fixed. In contrast, a radius variation of 
the tubular membrane accompanies changes in the tube 
volume and the area difference between the two leaflets. 
Therefore, it has to take the Laplace pressure into ac- 
count or requires an additional numerical technique to 
exchange the solvent particles or lipids between the up- 
per and lower sides of the bilayers. The new buckling 
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FIG. 3: (Color online) Area Axy dependence of the surface 
tension for Kcv/fceT = 9, 21, or 44 and = 800. The length 
Lx is varied with constant Ly — 24(t (x,o, A) or constant 
ratio Ly — Lx/2 (o, □, V). The solid lines are given by Eqs. 
(dl and ((9|. The dashed lines are given by the approximation, 
Eqs. (fTO )) and (fTTj) . 



method is suitable for measuring the bending rigidity of 
these membranes. 

We neglect the effects of thermal fluctuations in our 
analysis. The excess area induced by the thermal undu- 
lations shows a slight increase with deceasing area in the 
buckling simulation. However, this does not have a large 
effect on the surface tensions. For a squared membrane 
Lx = Ly, the thermal fluctuations can induce a flip be- 
tween the buckling in the x and y directions slightly be- 
low the buckling transition point. In the meshless mem- 
brane at Lx = Ly and N = 1600, this flip is observed for 
Axy/Na'^ > 1.3. Figure 2] shows the time development of 
the aspect ratio of the buckled amplitude. 



\h{qxi)\' + \h[qyi)[^ 



(12) 



where h{qxi) = z.j eyiY>{-2TTixj / L^) and h{qyi) = 
Zj exTp{—2TTiyj / Ly). The membrane changes the 
buckling direction at A^y/Na"^ — 1.3, while not at 
Axy/Na'^ = 1.2. 
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FIG. 4: (Color online) Time development of the aspect ratio 
/laap of the buckled amplitude for A^y/Na^ = 1.2 and 1.3 at 
Kcv/feT = 21, Lx = Ly, and iV = 1600. 




FIG. 5: (Color online) Snapshots of multiple buckled mem- 
branes at A' = 6400, A^mb = 8, Ly/a = 24, Kcv/fer = 21, 
and dn-i/u ~ 10. Left and right panels show the mem- 
branes with complete and partial phase synchronization at 
AxyNmh/Na^ — 1.1 and 1.25, respectively. 



IV. MULTIPLE BUCKLED MEMBRANES 

The membrane thermal undulations generate entropic 
repulsive force / oc d^^ between tensionless fluid mem- 
branes with neighboring membrane distance d [ssj. Here, 
we consider the interactions between the buckled mem- 
branes. Since the buckled membranes have greater height 
amplitudes than the tensionless membrane, the mem- 
branes have stronger short-range repulsion with decreas- 
ing area A^y (see Fig. [S|). Then, the fluctuation ampli- 
tude of the neighboring membrane distance ((d/dm — 1)^) 
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decreases for the fixed mean distance dm = Lz/N^^^ 
where A^mb is the number of the membranes (see the 
dashed line in Fig. EJa)). 

Along with this reduction in the distance fluctua- 
tion, it is found that the buckling of the neighbor- 
ing membrane becomes synchronized in phase. The 
phase (j) is calculated from the lowest Fourier mode 
of the membrane height, h[qxi) — /lam exp(— — 
^■ZjexTp{—2Trixj/Lx), for each membrane. The phase 
difference 0nb between neighboring membranes ap- 
proaches zero as A^y decreases. The phase deviation of 
the neighboring membrane (?!'nb)/'/'ran is shown as a sohd 
line in Fig. IH^a), where (j)'^^ is normalized by the average 
for the random distribution, = 0?an = tt^/S. Thus, 

the translational order of the buckled shape appears in 
the z direction. This ordering is not a discrete transition 
but a gradual change, since it is a quasi-one-dimensional 
system. As A^y decreases, clusters of the synchronized 
membrane appear, and then all of membranes are syn- 
chronized at AxyNynh/Na-^ < 1.1. (see the snapshots in 
Fig. [5] and movies in EPAPS [111)- The interactions be- 
tween the membranes little change their surface tension 
in the simulated area range. As the mean distance dm 
increases, the synchronization requires smaller A^y (see 
Fig. mjb)). This synchronized buckling may act as a 
nucleus to form multi-lamellar vesicles in shear flow. 



V. SUMMARY 

We studied the elastica of fluid membrane. The buck- 
led shape and surface tension parallel to the buckling (x) 
direction are expressed by the formula used for the elas- 
tica of solid sheets. However, unlike the solid sheets, the 
surface tension of the fluid membranes in the perpen- 
dicular (y) direction shows large increases for decreasing 
projected area A^y. Additionally, multi-lamellar buckled 
membranes were found to have phase synchronization. 

The anisotropic surface tension would also appear for 
the buckled Langmuir monolayers in the fluid phase. 
It can be experimentally checked, if one separately 
measures the stress in two lateral directions. These 
anisotropy and synchronization may play a role in the 
breakup of the membrane under external fields. Recent 
experiments show that collagen-containing tubular vesi- 
cles exhibit elastica-shape under magnetic field [3^1 • It 
will be interesting to investigate the coupling of mechan- 
ical and external-field induced bucklings. 
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FIG. 6: (Color online) Phase synchronization of buckled mem- 
branes at iV = 6400, iVmb = 8, Ly/a = 24, and Kcv/fcBT = 21. 
(a) The area A^y dependence of the fluctuation amplitudes of 
the distance {{d/dm ~ 1)^) (□) and the phase difi'erence (<^nb) 
(o) between neighboring membranes at dm/cr = 10. (b) The 
phase difference (^nb) dependence on the mean distance dm 
between neighboring membranes at A^yNmh/Na^ — 1.1 (A) 
and 1.2 (o). 



Appendix A: Details of simulation model and 
method 



A membrane consists of N particles, which possess no 
internal degrees of freedom. The particles interact with 
each other via a potential 



C/ = £(C/rcp + C/att)+A:aC/a, 



(Al) 
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which consists of a repulsive soft-core potential U^cp with 
a diameter cr, an attractive potential J7att, and a curva- 
ture potential All three potentials only depend on 
the positions of the particles. In this paper, we employ 
the curvature potential based on the first-order moving 
least-squares (MLS) method (model II in Ref. [lO). We 
briefiy outline here the simulation technique, since the 
membrane model is explained in more detail in Ref. [20l 
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Curvature Potential 



with B = 0.126, and a C°°-cutofr function m 



A Gaussian function with C° 
as a weight function 



cutoff ^ is employed 







t) in,j < rc 



(A2) 



where r^j is the distance between particles i and j. This 



function is smoothly cut off at 
the parameters n = 12, = 3ct, and 

The degree of deviation from a plane, 
defined as 



ga 



9Dw 



. We use here 
= 1.5cr. 

"aplanarity" is 



(A3) 



9Ai A2A3 



(Ai + A2 + A3)(AiA2 + A2A3 + A3A1) ' 

where Ai, A2, and A3 are the eigenvalues of the 

weighted gyration tensor, aa/j — J^ji'^j ~ '^G)il3j — 

/3G)wmis(rij), where a,/3 = x,y,z and rc = 

Yl,iT^jWmis{rij)/ J2j'Wm\sin,j). The aplanarity can be 

calculated from three invariants of the tensor: and 

are determinant and trace, respectively, and is 

the sum of its three minors, M„ = a^xO-yy + o^yyCizz + 
222 

0-zzO,xx ~ 0,xy ~ '^yz ~ ^zx' 

The aplanarity ttpi takes values in the interval [0, 1] 
and represents the degree of deviation from a plane. 
This quantity acts like Ai for Ai ^ A2,A3, since api ~ 
9Ai/(A2 + A3) in this limit. Therefore, we employ the 
curvature potential 



(A4) 



where Q!pi(r,i) = when the i-th particle has two or less 
particles within the cutoff distance r. 



This po- 
tential increases with increasing deviation of the shape 
of the neighborhood of a particle from a plane, and fa- 
vors the formation of quasi-two-dimensional membrane 
aggregates. 



/cut(s) 



r exp{A(l 



(|s|/Scut)" 



-)} (S < Scut) 



(A6) 



1 (s > Scut) 

is employed. All orders of derivatives of /cut(s) are con- 
tinuous at the cutoff. In Eq. (IA5p . we use the parameters 
n = 12, j4 = 1, and Scut — 1-2. 

A solvent-free membrane requires an attractive inter- 
action which mimics the "hydrophobic" interaction. We 
employ a potential J/att of the local density of particles. 



Pi ^^fcutirij/a), 



(A7) 



with the parameters n — 12, Shaif = 1-8, and Scut = 
2.1. The factor A in Eq. (|A7p is determined such that 
/cut(shaif) = 0.5, which implies A = ln(2){(scut/shaif)" - 
1}. Here, pi is the number of particles in a sphere whose 
radius is approximately ratt = ShaifC. The potential U^tt 
is given by 

C/att = 0.25 ln[l + exp{-4(p, - p*)}] ~ C, (A8) 



where C = 0.25 ln{l + exp(4/9*)}. For pi < p* , the po- 
tential is approximately [/att — —Pi and therefore acts 
like a pair potential with C/att - -J2t<j'^fcutin,j/(T). 
For Pi > p*, this function saturates to the constant — C. 
Thus, it is a pairwise potential with cutoff at densities 
higher than pi > p*. In this paper, we use e/ksT = 4 
and p* = 6. 



3. Dynamics 

The buckling of the membrane is simulated by Brow- 
nian dynamics (molecular dynamics with Langevin ther- 
mostat). The motion of particles is determined by the 
underdamped Langevin equations 



2. Attractive and Repulsive Potentials 



dU ^dxi_ 
dTi dt 



(A9) 



The particles interact with each other in the quasi- two- 
dimensional membrane surface via the potentials C/rcp 
and f/att- The particles have an excluded- volume inter- 
action via the repulsive potential 

C/rcp = J2 cxp(-20(r,j/a - 1) + B)/cut(r-,;,,/a), (A5) 

i<j 



where m is the mass of a particle and C the friction 
constant. gi{t) is a Gaussian white noise which obeys 
the fluctuation-dissipation theorem. We employ the time 
unit T = C^a^ /k^T with m — ^r. The Lan gevi n equations 
are integrated by the leapfrog algorithm ]36l. Wf\ with a 
time step of M = 0.005r. 
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